% deflation probabilities
clear

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
year = 1850:1:2013;
%select_years = [1876, 1909, 1921, 1980, 1998, 2009];
select_years = [27, 60, 72, 130, 149, 160];
NY = size(select_years,2);
HZN = 10; % longest forecast horizon

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% catalog of data files
DFILE(1,:) = ['swuc_wpi_ar1_01'];
DFILE(2,:) = ['swuc_wpi_ar1_02'];
DFILE(3,:) = ['swuc_wpi_ar1_03'];
DFILE(4,:) = ['swuc_wpi_ar1_04'];
DFILE(5,:) = ['swuc_wpi_ar1_05'];
DFILE(6,:) = ['swuc_wpi_ar1_06'];
DFILE(7,:) = ['swuc_wpi_ar1_07'];
DFILE(8,:) = ['swuc_wpi_ar1_08'];
DFILE(9,:) = ['swuc_wpi_ar1_09'];
DFILE(10,:) = ['swuc_wpi_ar1_10'];
DFILE(11,:) = ['swuc_wpi_ar1_11'];
DFILE(12,:) = ['swuc_wpi_ar1_12'];
DFILE(13,:) = ['swuc_wpi_ar1_13'];
DFILE(14,:) = ['swuc_wpi_ar1_14'];
DFILE(15,:) = ['swuc_wpi_ar1_15'];
DFILE(16,:) = ['swuc_wpi_ar1_16'];
DFILE(17,:) = ['swuc_wpi_ar1_17'];
DFILE(18,:) = ['swuc_wpi_ar1_18'];
DFILE(19,:) = ['swuc_wpi_ar1_19'];
DFILE(20,:) = ['swuc_wpi_ar1_20'];
DFILE(21,:) = ['swuc_wpi_ar1_21'];
DFILE(22,:) = ['swuc_wpi_ar1_22'];
DFILE(23,:) = ['swuc_wpi_ar1_23'];
DFILE(24,:) = ['swuc_wpi_ar1_24'];
DFILE(25,:) = ['swuc_wpi_ar1_25'];
DFILE(26,:) = ['swuc_wpi_ar1_26'];
DFILE(27,:) = ['swuc_wpi_ar1_27'];
DFILE(28,:) = ['swuc_wpi_ar1_28'];
DFILE(29,:) = ['swuc_wpi_ar1_29'];
DFILE(30,:) = ['swuc_wpi_ar1_30'];
DFILE(31,:) = ['swuc_wpi_ar1_31'];
DFILE(32,:) = ['swuc_wpi_ar1_32'];
DFILE(33,:) = ['swuc_wpi_ar1_33'];
DFILE(34,:) = ['swuc_wpi_ar1_34'];
DFILE(35,:) = ['swuc_wpi_ar1_35'];
DFILE(36,:) = ['swuc_wpi_ar1_36'];
DFILE(37,:) = ['swuc_wpi_ar1_37'];
DFILE(38,:) = ['swuc_wpi_ar1_38'];
DFILE(39,:) = ['swuc_wpi_ar1_39'];
DFILE(40,:) = ['swuc_wpi_ar1_40'];

NF = size(DFILE,1);
NB = 21;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  compute RA and QA
load(DFILE(NB,:),'vD')
[N,T,P] = size(vD);
RA = zeros(T,(NF-NB+1)*P);
QA = zeros(T,(NF-NB+1)*P);
RA(:,1:P) = vD(1,:,:);
QA(:,1:P) = vD(2,:,:);
clear vD 
for i = NB+1:NF,
  j = i - NB + 1; 
  load(DFILE(i,:),'vD')
  RA(:,(j-1)*P+1:j*P) = vD(1,:,:);
  QA(:,(j-1)*P+1:j*P) = vD(2,:,:);
  clear vD;
end
NMC = size(RA,2); % draws in posterior

% load WA
load(DFILE(NB,:),'WD')
[P,N,N] = size(WD);
WA = zeros((NF-NB+1)*P,N,N);
WA(1:P,:,:) = WD;
clear WD 
for i = NB+1:NF,
  j = i - NB + 1; 
  load(DFILE(i,:),'WD')
  WA((j-1)*P+1:j*P,:,:) = WD;
  clear WD;
end

% compute SA
load(DFILE(NB,:),'SD')
[P,K,T] = size(SD);
SA2 = zeros((NF-NB+1)*P,T); % mu
SA2(1:P,:) = SD(:,2,:);
clear SD 
for i = NB+1:NF,
  j = i - NB + 1; 
  load(DFILE(i,:),'SD')
  SA2((j-1)*P+1:j*P,:) = SD(:,2,:);
  clear SD;
end

% random subsample from the MCMC sample
ords = sort(round(NMC*rand(10000,1)));
NMC2 = size(ords,1);
SA2s = SA2(ords,:); clear SA2;
WAs = WA(ords,:,:); clear WA;
RAs = RA(2:T+1,ords); clear RA;
QAs = QA(2:T+1,ords); clear sA;

deflation_prob = zeros(HZN,NMC2);
DPR = zeros(2,NMC2,NY);

IQR = zeros(2,3,T);
matlabpool local 8
for yy = 1:T, 
    yy 
    parfor ii = 1:NMC2, 
        deflation_prob(:,ii) = posterior_predictive_distribution(squeeze(WAs(ii,:,:)),RAs(yy,ii),QAs(yy,ii),SA2s(ii,yy),HZN,10000);
    end
    ss = sort(deflation_prob([5 10],:)');
    IQR(:,:,yy) = ss([.25 .5 .75]*NMC2,:)';
    
    % store deflation probabilities for peak and trough years
    [mm,ind] = min(abs(yy*ones(1,NY) - select_years));
    if mm == 0,
        DPR(:,:,ind) = deflation_prob([5 10],:);
    end
end
matlabpool close

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% median and iqr for deflation probabilities
figure
subplot(2,2,1)
plot(year,squeeze(IQR(1,2,:)),'-k','Linewidth',2); hold on;
plot(year,squeeze(IQR(1,[1 3],:)),':k','Linewidth',2); hold off;
title('5 years ahead','Fontsize',16)
legend('Median','Interquartile range','Location','NorthEast')
ylabel('Prob( p_{t+h} - p_{t} < 0)','Fontsize',16)
axis([-inf inf 0 1])
set(gca,'Fontsize',16)
grid

subplot(2,2,2)
plot(year,squeeze(IQR(2,2,:)),'-k','Linewidth',2); hold on;
plot(year,squeeze(IQR(2,[1 3],:)),':k','Linewidth',2); hold off;
title('10 years ahead','Fontsize',16)
axis([-inf inf 0 1])
set(gca,'Fontsize',16)
grid


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% comparing post WWII with pre WWI: deflation probabilities
EDPR = squeeze(mean(DPR,2)); % mean smoothed DPR
rm5 = EDPR(1,:)/EDPR(1,1); % 1876 is base year
rm10 = EDPR(2,:)/EDPR(2,1);
[rm5(4:6)' rm10(4:6)']
% figure
% subplot(2,3,1)
% bar([rm5(4:7)' rm10(4:7)'],0.3); hold on;
% plot(ones(size(rm5(4:7)')),'--k','Linewidth',2); hold off
% title('Relative to 1876','Fontsize',14)
% legend('5 years ahead','10 years ahead','Location','Northwest')
% ylabel('Ratio of Mean DPR','Fontsize',14)
% axis([-inf inf 0 1.2])
% set(gca,'Fontsize',14)
% grid
% 
rm5 = EDPR(1,:)/EDPR(1,2); % 1906 is base year
rm10 = EDPR(2,:)/EDPR(2,2);
[rm5(4:6)' rm10(4:6)']
% subplot(2,3,2)
% bar([rm5(4:7)' rm10(4:7)'],0.3); hold on
% plot(ones(size(rm5(4:7)')),'--k','Linewidth',2); hold off
% title('Relative to 1907','Fontsize',14)
% axis([-inf inf 0 1.2])
% set(gca,'Fontsize',14)
% grid

rm5 = EDPR(1,:)/EDPR(1,3); % 1921 is base year
rm10 = EDPR(2,:)/EDPR(2,3);
[rm5(4:6)' rm10(4:6)']
% subplot(2,3,3)
% bar([rm5(4:7)' rm10(4:7)'],0.3); hold on;
% plot(ones(size(rm5(4:7)')),'--k','Linewidth',2); hold off
% title('Relative to 1931','Fontsize',14)
% axis([-inf inf 0 1.2])
% set(gca,'Fontsize',14)
% grid

% PROBABILITY OF A DECLINE RELATIVE TO A BASE YEAR
CDPR5 = squeeze(DPR(1,:,:));
CDPR10 = squeeze(DPR(2,:,:));
PR = zeros(3,2);
for yy = 4:6,
    rcdpr = CDPR5(:,yy)./CDPR5(:,1); 
    DD5 = zeros(size(rcdpr));
    DD5(rcdpr<1) = 1;
    PR(yy-3,1) = mean(DD5');

    rcdpr = CDPR10(:,yy)./CDPR10(:,1);
    DD10 = zeros(size(rcdpr));
    DD10(rcdpr<1) = 1;
    PR(yy-3,2) = mean(DD10');
end
PR
% subplot(2,3,4)
% bar(PR,0.3); hold on;
% plot(.9*ones(size(PR)),'--k','Linewidth',2); hold on;
% plot(.1*ones(size(PR)),'--k','Linewidth',2); hold off;
% ylabel('Probability of decline','Fontsize',14)
% axis([-inf inf 0 1.01])
% set(gca,'Fontsize',14)
% grid

PR = zeros(3,2);
for yy = 4:6,
    rcdpr = CDPR5(:,yy)./CDPR5(:,2); 
    DD5 = zeros(size(rcdpr));
    DD5(rcdpr<1) = 1;
    PR(yy-3,1) = mean(DD5');

    rcdpr = CDPR10(:,yy)./CDPR10(:,2);
    DD10 = zeros(size(rcdpr));
    DD10(rcdpr<1) = 1;
    PR(yy-3,2) = mean(DD10');
end
PR
% subplot(2,3,5)
% bar(PR,0.3); hold on;
% plot(.9*ones(size(PR)),'--k','Linewidth',2); hold on;
% plot(.1*ones(size(PR)),'--k','Linewidth',2); hold off;
% axis([-inf inf 0 1.01])
% set(gca,'Fontsize',14)
% grid

PR = zeros(3,2);
for yy = 4:6,
    rcdpr = CDPR5(:,yy)./CDPR5(:,3); 
    DD5 = zeros(size(rcdpr));
    DD5(rcdpr<1) = 1;
    PR(yy-3,1) = mean(DD5');

    rcdpr = CDPR10(:,yy)./CDPR10(:,3);
    DD10 = zeros(size(rcdpr));
    DD10(rcdpr<1) = 1;
    PR(yy-3,2) = mean(DD10');
end
PR
% subplot(2,3,6)
% bar(PR,0.3); hold on;
% plot(.9*ones(size(PR)),'--k','Linewidth',2); hold on;
% plot(.1*ones(size(PR)),'--k','Linewidth',2); hold off;
% axis([-inf inf 0 1.01])
% set(gca,'Fontsize',14)
% grid
